function [tspeed, tvel_x, tvel_y, tshrink]=xyt2speed(tx,ty,tt,ta,filt)
%
% [tspeed tvel_x tvel_y tshrink]=xyt2speed(tx,ty,tt,[ta],[filt])
% tx = x position
% ty = y position
% tt = time vector (monotonically increasing is best)
% ta = area [optional]
% all must be vectors and of the same length.
%
% filt = filtfilt kernel size for box filter OR send FiltObj.Numerator
% FiltSpec = fdesign.lowpass('N,Fc',50,50,500); 50 point, 50Hz cutoff filter for 500hz samplerate data
% FiltObj = design(FiltSpec);
%
% coe@queensu.ca
if nargin < 3
    help(mfilename)
    return
end
if nargin < 5
    filt=0;
end
if nargin < 4
    ta=[];
    tshrink=[];
end
if length(ta)~=length(tt)% easter egg to support three inputs ignoring area
    filt=ta;
    ta=[];
    tshrink=[];
end
filt=double(filt);
if length(filt)==1 % to support different kernel types
    if filt==1
        filt=2;
    end
    filt=ones(filt,1)/filt;% make a box kernel here
end
if filt
    if sum(single(filt))~=1
        warning('xyt2speed Filter SUM =%f',sum(filt)) %#ok<WNTAG>
    end
end
tvel_x=single(nan(length(tx),1));
tvel_y=tvel_x;
gooddata=~isnan(tx)&~isnan(ty);
if sum(gooddata)<length(tt)/2
    tspeed=zeros(size(tt));
    tvel_x=tspeed;
    tvel_y=tspeed;
    tshrink=tspeed;
    warning(' eye data is too poor to use in xyt2speed') %#ok<WNTAG>
    return
end
tx2=tx(gooddata);
ty2=ty(gooddata);
tt2=tt(gooddata);
EE=length(tt2);
tdiff=[(tt2(2)-tt2(1));(tt2(3:EE)-tt2(1:EE-2))/2;(tt2(EE)-tt2(EE-1))];
xdiff=[(tx2(2)-tx2(1));(tx2(3:EE)-tx2(1:EE-2))/2;(tx2(EE)-tx2(EE-1))];
ydiff=[(ty2(2)-ty2(1));(ty2(3:EE)-ty2(1:EE-2))/2;(ty2(EE)-ty2(EE-1))];
tvel_x(gooddata)=xdiff./tdiff*1000; % deg per sec
tvel_y(gooddata)=ydiff./tdiff*1000; % deg per sec
if filt
    tvel_x(gooddata)=single(filtfilt(filt,1,double(tvel_x(gooddata))));
    tvel_y(gooddata)=single(filtfilt(filt,1,double(tvel_y(gooddata))));
end
tspeed=sqrt(tvel_x.^2+tvel_y.^2);
if ~isempty(ta) % nargin==4 && nargout==4
    tshrink=zeros(length(tx),1);
    ta2=ta(gooddata);
    adiff=[(ta2(2)-ta2(1));(ta2(3:EE)-ta2(1:EE-2))/2;(ta2(EE)-ta2(EE-1))];
    if filt
        adiff=single(filtfilt(filt,1,double(adiff)));
    end
    tshrink(gooddata)=abs(adiff./tdiff);
end
return
